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ABSTRACT 

This  report  describes  the  main  features  of  a  digital 
computer  program  for  calculating  the  pressure  distribution  along 
bodies  in  two-dimensional  cavity  flow  in  an  infinite,  incompres¬ 
sible,  inviscid  fluid.  The  pressure  distribution  is  used  to 
compute  drag,  lift,  and  moment  coefficients.  The  program  is 
based  on  a  theory  of  cavitated  flows  developed  in  recent  years 
by  Wu.  In  the  present  report,  the  theory  is  applied  to  flow  con¬ 
figurations  characterized  by  two  fixed  detachment  points  and  to 
sections  whose  shape  in  the  wetted  region  is  arbitrary,  apart  from 
the  restriction  that  there  be  no  sharp  corners  interior  to  the 
wetted  region.  Calculations  are  presented  for  parabolic  sections 
of  10-to-l  fineness  ratio  and  1-to-l  fineness  ratio  at  several 
angles-of -attack. 


NOMENCLATURE 


A’ 

B 

B' 

C 

C* 


=  upper  detachment  point  in  actual  flow;  also,  a 

real  amplitude  associated  with  the  complex  velocity 
potential 

=  upper  detachment  point  in  basic  flow 

=  lower  detachment  point 

=  lower  detachment  point  in  basic  flow 

=  point  of  intersection  of  ^-axis  with  base  of  body 

-  point  of  intersection  of  'x-axis  with  base  of  wedge 

in  basic  flow 

P-PC 


*  pressure  coefficient  (normalized  pressure)  = 

F. 

=  drag  coefficient  =*  * — =_ 

7  pU  L 


PS‘PC 


-  lift  coefficient  = 


c 


m 


g(t,tQ) 

h(w) 


I(n)(t') 

i 


moment  coefficient 


M 

T  ■"'^2 

^  pU  L 


component,  in  the  direction  of  the  incident  stream, 
of  the  hydrodynamic  force  per  unit  span 

component  of  same  force  normal  to  the  incident  stream 

form  of  complex  velocity  potential  in  the  t-plane 

X 

functional  relationship  between  y-coordinate  of  point 
wetted  perimeter  and  auxiliary  variable,  w 

arc-length  integral  in  w,  extending  from  w^  to 

integrand  of  t* integral  in  Equation  (3.5) 
subscript  denoting  iC^  mqsh  point;  i  *  1,2,  .... 


vii. 


j(n)(t')  _  x-integral  In  Equation  (3.5) 


N(t) 


(n>- 


=  subscript  denoting  j  mesh  or  Simpson  point; 
j  =  1,2,  .... 

=  reference  length  (=*  &C  of  Figure  2.1)  in  units  of 
which  all  normalized  lengths  are  expressed 

=»  length  of  side  of  wedge  in  basic  flow 

=*  hydrodynamic  moment  about  point  0 

*  j  s(t) 

=  number  of  fine -mesh  intervals  in  w  or  t 

=  number  of  rough-mesh  intervals  in  w  or  t  on  upper 
profile 

=  number  of  rough -mesh  intervals  in  t  on  lower  profile 
in  the  n  iteration 

=  total  number  of  mesh  points  in  t  in  nc^  iteration 

=«  superscript  denoting  n^  iteration;  n*0,l,  .... 

a  nose  of  body,  and  origin  of  coordinate  system 

»  mesh  point  on  upper  profile  separating  rough  from 
fine  mesh 

»  provisional  point  on  lower  profile  separating  fine 
from  rough  mesh 

<«  final  mesh  point  on  lower  profile  separating  fine 
from  rough  mesh 

«  pressure  anywhere  on  wetted  perimeter 
a  cavity  pressure 

*  stagnation  pressure 
«  pressure  at  infinity 

»  real  number  such  that  0  £  r  <  1 

*  total  length  of  wetted  perimeter 

$ 

"  i 


viii. 


s(t) 

3(n)(t) 

s  <»-D 
8j 

s[w] 


arc  length  measured  from  upper  detachment  point  A 
to  any  point  on  wetted  perimeter 


$ 

L 


s  as  a  function  of  t 

approximation  to  s(t)  obtained  after  n^  iteration 


8<n'1>(Tj,n> 

s  as  a  function  of  w 


stWjjj.] 

s[wFg) 

variable  -  in  general  complex  -  into  whose  plane 
the  physical  plane  is  mapped  so  that  the  complex 
velocity  potential  takes  the  form  shown  in  Equation 
(2.1).  However,  in  all  functions  and  Integrals  in 
this  paper,  t  takes  on  only  real  values  and  is  con¬ 
fined  to  the  interval  -1  £  t  £  1. 

image  in  the  t-plane  of  the  point  at  infinity  in 
the  physical  plane 

approximation  to  tQ  obtained  after  n^  iteration, 
n  ®  0,1,  .... 

i  mesh  point  in  t  in  n  iteration,  i  *  1,2, 
n  "**  0,1,  .... 

value  of  t  satisfying  equation,  s^n_1^(t)  -  s^ 

value  of  t  satisfying  equation,  s^n“^(t)  *  spR 

the  particular  t^n^  lying  closest  to  tpR^* 
mesh  width  in  t  in  basic  flow 

width  of  rough  mesh  corresponding  to  the  n  iteration 
width  of  fine  mesh  corresponding  to  the  n  iteration 


ix 


U 

V 
w 

wi 

w] 

w, 

w, 

w 

1 

wmax 

<AU>R 

(Aw)f 

$ 


speed  of  incident  stream 


RF 


FR 


FR 


min 


RF 

★ 

*FR 

*FR 

♦ 

$ 


e 


-  (l-NT) 


1 


=  ±vr=  auxiliary  vrriable  used  in  representing  the 
wetted  profile 

=  the  ifck  mesh  point  in  w,  i  =  1,2,  .... 

:  v#7 

VXFR 

=  the  particular  w^  lying  closest  to  w^R 
*  W1 

=  Wnr+nf+nr+1 

s»  width  of  rough  mesh  in  w 
=*  width  of  fine  mesh  in  w 


abscissa  of  arbitrary  point  on  wetted  perimeter 
x 

i 


=*  x-coordinate  of  P, 


RF 


x-coordinate  of  P, 


FR 


*RF 

T 

^FR 

nr- 

L 


TR 


T(») 


ordinate  of  arbitrary  point  on  wetted  perimeter 

$  +  i$ 

angle  between  Incident  flew  direction  and  positive 
$-axis 

angle,  measured  positive  counterclockwise,  between 
tangent  vector  to  wetted  perimeter  and  positive 
$-axis.  Tangent  vector  pointy  in  direction  of  in¬ 
creasing  s. 

3  expressed  as  a  function  of  s 


X. 


relative  error 

fraction  of  v  by  which  the  half-angle  of  the  basic 
flow  wedge  is  measured 

density  of  the  fluid 

.  .  PC"Poo 

cavitation  number  =  -  — — r~ 

Ps_p, 

real  integration  variable;  -1  £  t  £  1 

an'  Simpson  point  associated  with  the  x-integral 
of  Equation  (35)  Me  n  iteration. 
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1.  INTRODUCTION  AND  SUMMARY 

As  part  of  an  investigation  into  ventilated  flow  past 
surface-piercing  struts,  TRG  has  been  developing  a  digital  com¬ 
puter  program  for  calculating  the  pressure  distribution  along 
bodies  in  two  dimensional  cavity  flow  in  an  infinite,  incompres¬ 
sible,  inviscid  fluid.  From  this  pressure  distribution  we  compute 
the  drag  coefficient  and  -  in  the  case  of  non-zero  angle  of  attack- 
the  lateral  or  lifting  force  coefficient  and  the  hydrodynamic 
moment  coefficient.  The  program  is  based  on  a  theory  of  such 
flovc  developed  in  recent  years  by  As  of  this  writing, 

we  nave  applied  the  theory  to  flow  configurations  characterized 
by  two  fixed  detachment  points  and  to  sections  whose  shape  in  the 
wetted  region  is  arbitrary,  apart  from  the  restriction  that  there 
be  no  sharp  corners  interior  to  the  wetted  region.  To  the  best  of 
our  knowledge,  Wu' s  theory  has  not  previously  been  applied  to  such 
problems.  Future  work  may  include  situations  in  which  one  or  both 
detachment  points  are  smooth  and/or  sharp  corners  interior  to  the 
wetted  region  are  permitted. 

In  this  report  we  describe  the  main  features  of  our  pro¬ 
gram,  and  we  exhibit  and  discuss  the  computed  results  obtained 
thus  far  in  a  number  of  cases.  It  will  be  seen  that  in  two  of 
these  cases  our  computational  scheme  worked  well,  in  one  case  it 
gave  results  whose  interpretation  we  have  not  yet  decided  upon,  and 
in  one  case  it  was  unable  to  yield  meaningful  answers.  The  prin¬ 
ciple  cause  of  difficulty  was  probably  the  insufficient  accuracy 
attained  in  the  evaluation  of  certain  integrals;  any  further  work 
on  the  two  dimensional  problem  should  therefore  include  an  effort 


to  improve  this  accuracy.  There  is  no  obstacle  other  than  pos¬ 
sibly  increased  running  time,  to  achieving  this  higher  accuracy. 


2. 


Before  going  into  details,  we  summarize  the  results. 
Allusions  to  convergence  refer  to  an  iteration  procedure  which  is 
one  of  the  principal  features  of  the  computational  scheme,  and  which 
is  described  in  Section  3. 

Case  1.  A  10-to-l  fineness  ratio  parabola  at  zero  angle- 
of-attack  and  zero  cavitation  number.  (See  Figure  1.1.)  The  se¬ 
quence  of  computed  pressure  distributions  does  net  converge.  This 
may  be  seen  from  the  fact  that  the  first  three  pressure  distribu¬ 
tions  in  the  sequence  are  almost  symmetric  with  respect  to  the  nose 
of  the  parabola  (the  true  distribution  is  symmetric)  while  the 
fourth  distribution  is  close  to  antisymmetric.  (See  Figure  4.1.) 

The  corresponding  values  of  drag  coefficient  fluctuate  from  one 
iteration  to  the  next  in  what  seems  to  be  a  random  manner,  with  a 
maximum  deviation  of  about  20%  above  and  below  the  mean  for  the 
four  iterations. 

Case  2.  A  1-to-l  parabola  at  zero  angle-of-attack  and 
zero  cavitation  number.  (See  Figure  1.2.)  Evidence  that  the  se¬ 
quence  of  computed  pressure  distributions  converges  is  provided 
by  the  distributions  arising  from  the  sixth  and  seventh  iterations, 
which  differ  by  no  more  than  about  3.5  percent  anywhere  on  the 
wetted  surface,  while  the  corresponding  drag  coefficients  agree  to 
better  than  a  small  fraction  of  a  percent.  Confirmation  of  the 
correctness  of  the  final  (seventh)  distribution  is  provided  by  a 
comparison  with  a  simple  approximate  algebraic  expression  for  the 

f  31 

pressure  along  a  parabola,  derived  by  Johnson1,  J.  The  agreement 


CAVITY 


Figure  1.1 

10-to-l  PARABOLA,  0°  ANGLE -OF -ATTACK 


Figure  1.2 

1-to-l  PARABOLA,  0°  ANGLE -OF -ATTACK 
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between  our  result  and  compucations  based  on  Johnson’ s  formula  is 
quite  close. 

Case  3.  A  1-to-l  parabola  at  1°  angle-of-attack  and  zero 
cavitation  number.  (See  Figure  1.3.)  Convergence  is  suggested  by 
the  fact  that  the  difference  between  successive  pressure  distribu¬ 
tions  diminishes  with  iteration  number  over  the  entire  wetted  sur¬ 
face.  However,  the  rate  of  convergence  is  slow,  and  even  the  thir¬ 
teenth  and  fourteenth  distributions  differ  by  as  mucn  as  15  percent 
at  a  point  near  the  nose  where  the  pressure  is  40  percent  of  its 
stagnation  value.  Further  away  from  the  nose,  where  the  pressure 
is  lower,  the  difference  in  the  distributions  rises  to  25  percent. 
The  drag  coefficient  remains  constant  over  all  fourteen  iterations, 
apart  from  a  random  fluctuation  whose  rms  value  is  less  than  1  per¬ 
cent.  The  fluctuation  is  probably  caused  by  numerical  inaccuracies. 
These  inaccuracies  have  a  more  serious  effect  when  computing  the 
lateral  force  and  moment  coefficients,  in  which  they  produce  consider 
ably  larger  fluctuations.  The  reason  for  this  will  be  discussed  in 
Section  4. 

Case  4.  A  1-to-l  parabola  at  a  5°  angle-of-attack  and 
zero  cavitation  number.  (See  Figure  1.4w)  In  this  case,  after 
fourteen  iterations  there  is  no  indication  of  convergence  to  a  limit 
ing  pressure  distribution.  On  the  other  hand,  all  of  the  distribu¬ 
tions  are  similar  in  form.  Moreover,  apart  from  the  third  distri¬ 
bution,  which  has  a  large  negative  dip,  the  pressure  distributions 
fall  within  a  range  of  values  such  that  any  one  of  them  appears 
physically  reasonable.  This  suggests  that  the  failure  to  converge 
is  a  numerical  phenomenon,  an  idea  that  is  reinforced  by  the  fact 


Figure  1.3 

1-to-l  PARABOLA,  1°  ANGLE -OF -ATTACK 


Figure  1.4 

1-to-l  PARABOLA,  5°  ANGLE- OF- ATTACK 
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that,  as  in  the  1°  case,  the  drag  coefficient  varies  in  what  ap¬ 
pears  to  be  a  random  manner  as  a  function  of  iteration  number.  Al¬ 
though  the  fluctuation  of  the  drag  coefficient  about:  its  mean  value 
is  larger  than  in  the  1°  case,  it  is  still  small,  having  an  rms 
value  of  4  to  5  percent.  Also  as  in  the  1°  case,  the  fluctuations 
of  the  lift  and  moment  coefficients  are  considerably  larger  than  those 
of  the  drag  coefficient,  for  reasons  to  be  discussed  in  Section  4. 

It  should  be  remarked  that  the  computation  time  in  all 
four  cases  averaged  to  somewhat  less  than  0.01  hours  per  iteration 
on  the  IBM  7094  computer. 

To  render  intelligible  the  computational  procedure  for 
calculating  the  pressure  distribution  on  the  body,  and  the  corre¬ 
sponding  drag,  lift,  and  moment  coefficients,  it  is  necessary  to 
give  a  brief  account  of  the  underlying  theory.  Accordingly,  the 
next  section  contains  such  an  account.  In  Section  3,  we  describe 
the  computational  scheme  itself.  Although  we  do  not  exhibit  the 
actual  computer  program,  we  do  explain  all  procedures  colled  for 
in  the  program  specification,  other  than  standard  subroutines.  In 
Section  4,  we  present  some  of  the  numerical  results  obtained  for 
the  aforementioned  cases.  Finally,  Section  5  lists  some  recommen¬ 
dations  for  future  improvements  in  the  computational  procedure. 
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2.  OUTLINE  OF  THEORY 

The  configuration  under  discussion  is  shown  in  Figure 

2.1. 

In  the  figure,  AOBA  represents  the  perimeter  of  the 
body  upon  which  is  incident  a  uniform  flow  of  complex  velocity 
Ue^a.  The  shaded  region  represents  the  cavity,  the  points  A  and 
B  being  the  fixed  detachment  points  from  which  emanate  the  two 
free  (constant  pressure)  streamlines  that  bound  the  cavity.  0, 
the  origin  of  the  coordinate  system,  is  generally  placed 
most  conveniently  at  the  point  of  maximum  curvature  of  the  peri¬ 
meter;  for  the  usual  bodies  this  point  is  well  defined.  It  is 
convenient,  also,  to  have  the  $-axis  tangent  to  the  perimeter  at 
0.  P(s)  is  any  running  point  on  the  perimeter,  where  §  is  the 
arc  length  from  A  to  P.  0  is  the  angle,  measured  positive  counter¬ 
clockwise,  between  the  positive  ^-direction  and  the  tangent  vec¬ 
tor  at  P.  This  vector  points  in  the  direction  of  increasing 

•jkf 

The  significance  of  the  points  P^,  Pj^,  and  PpR  with  abscissae 
^RF*  ^fr’  and  ^FR  explained  presently. 

As  mentioned,  our  computation  procedure  is  based  on  a 
theory  due  to  wherein  the  region  of  potential  flow 

exterior  to  the  body  and  cavity  in  the  physical  ($«3+iy)  plane 
is  mapped  conformally  into  the  interior  of  a  8 lifted  semicircle 
in  the  plane  of  a  new  complex  variable  denoted  by  t.  (Here,  t 
has  no  connection  with  the  time.)  This  mapping  has  the  property 

y  11  ~  — . 

In  Hhi' s  mapping,  the  origin  of  the  rectangular  coordinate  system 
in  the  physical  plane  is  at  A,  but  this  has  no  effect  on  the 
equations  needed  in  the  present  work. 


Figure  2.1 


GENERAL  CONFIGURATION  UNDER  DISCUSSION 


that  the  velocity  potential  in  the  complex  t  plane,  call  it 
f(t,tQ),  takes  the  form  of  an  elementary  algebraic  function  of 
t,  viz. , 

=  At2  £(t-t0)  (t-tQ)  (t-t"1)  (t-t'Sj  1  (2. 

where  tQ  is  the  image,  in  the  t -plane,  of  the  point  £=«>,  A  is 
a  real  constant  (not  to  be  confused  with  point  A  of  Figure  2.1) 
to  be  determined,  and  the  bar  (in  this  case)  denotes  the  complex 
conjugate. 

Naturally,  the  complex  function  that  accomplishes  the 
mapping  is  not  known  in  advance,  and  its  determination  consti¬ 
tutes  the  main  problem.  However,  if  as  in  our  case,  one  is  not 
attenpting  to  obtain  the  entire  flow  but  only  the  pressure  dis¬ 
tribution  on  the  body,  one  can  get  by  with  the  computation  of  a 
good  deal  less  than  the  entire  mapping  function. 

The  base  of  the  semi -circle  in  the  t -plane  is  the  in¬ 
terval  -1  £  t  £  1  on  the  real  t-axis,  and  this  interval  is  also 
the  image  in  the  t-plane  of  the  wetted  perimeter  AOB  of  the  body. 
Thus  as  t  goes  from  -1  to  1,  P  goes  from  A  to  B  via  0,  while  8 
varies  from  0  to  £,  where  $  is  the  length  of  the  wetted  perimeter. 
As  will  become  evident,  the  determination  of  the  unknown  function 
$(t)  is  an  essential  step  in  obtaining  the  pressure  distribution 
on  the  body. 

Before  writing  expressions  for  these  quantities,  we 
Introduce  certain  normalizations  into  the  problem.  Specifically, 
we  normalize  ail  lengths  with  respect  to  some  convenient  length, 
t.,  defined  in  the  plane  of  the  flow.  \  in  our  case  is  taken  to 
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r  %  a 

be  OC,  the  length  of  x-axis  intercepted  by  the  body.  Then,  if  s 

£ 

represents  normalized  arc  length,  s  =  Likewise,  all  other 
normalized  lengths  are  denoted  by  removing  the  hat  symbol  from 
the  corresponding  letter.  Further,  we  normalize  all  velocity 
magnitudes with  respect  to  the  constant  speed  of  the  fluid  along 
the  free  streamlines  that  bound  the  cavity.  Finally,  we  introduce 
three  dimensionless  coefficients,  viz,,  a  two-dimensional  drag 
coefficient,  c^,  a  two-dimensional  lateral  force  or  lift  coef¬ 
ficient,  c^,  and  a  two-dimensional  moment  coefficient,  cm.  These 
are  defined  as  follows: 


c 


Fd 

$5* 


(2.1*) 


*1 


(2.2b) 


(2'3) 

where  is  the  component,  in  the  direction  of  the  Incident 
stream,  of  the  hydrodynamic  force  per  unit  span  acting  on  the 
body,  F|  is  the  component  of  the  same  force  normal  to  the  stream, 
and  M  is  the  moment  (measured  positive  counterclockwise)  of  that 
force  about  the  nose,  0.  U  is  the  unnormalized  speed  of  the  in¬ 
cident  stream  and  p  is  fluid  density. 

Instead  of  the  pressure  itself,  it  is  convenient  to 
coelute  the  dimensionless  pressure  coefficient,  Cp,  defined  by 
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P-Pr 
C  =■• - 1 

p  ps'pc 


(2.4) 


where  p  is  the  fluid  pressure  anywhere  along  the  perimeter  of  the 
body,  pc  is  the  pressure  in  the  cavity,  and  pg  is  the  stagnation 
pressure.  It  can  be  shown  that  the  drag,  lift,  and  moment  co¬ 
efficients  are  then  expressible  in  terms  of  Cp  by  means  of  the 
relationships 


where 


c^  =  cx  cosa  +  Cy  sina 


c l  -  "  cx  sina  +  Cy  cosa 


cm  " 


(l+o)  ^  Cp(s)^x(s)cosF(s)+y(e)sinp(s)|  ds 


s 

l 


cx  -  -  (l+o)  I  Cp(s)  sin^(s)  ds 


/ 


Cy  »  (l+o)  J  Cp(s)  cos]5(s)  ds 


(2.5a) 

(2.5b) 

(2.5c) 

(26a) 

(2.6b) 


In  these  equations,  a  is  the  angle,  measured  positive 
counterclockwise,  between  the  incident  flow  velocity  and  the  posi¬ 
tive  x-axis,  S  is  the  normalized  length  of  the  wetted  perimeter 
and  a  is  the  cavitation  number,  defined  by 


pc’p- 

PS’P« 


(2.7) 


where  pw  is  the  fluid  pressure  at  infinity.  ]5(s)  denotes  the 
functional  dependence  of  the  angle  £  on  the  arc  length  s,  while 
x(s)  and  y(s)  denote  the  corresponding  dependences  of  the  x  and 
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y  coordinates  of  a  point  on  the  perimeter.  These  functional  de¬ 
pendences  are  known;  being  determined  solely  by  the  geometry  of 
the  body. 

Next  we  change  the  integration  variable  from  s  to  t. 

It  is  shown  in  Wu^’^J  that  Cp[t],  defined  as  Cp(s(t)),  can  be 
written 


In  Eq.  (2.8),  f(t,tQ)  is  the  image  in  the  t-plane  of 
the  normalized  velocity  potential  (see  Equation  (2.1)),  s(t)  is 
the  still-to-be-determined  normalized  arc  length  as  a  function  of 
t,  a;.d  the  subscript  t  denotes  differentiation  with  respect  to  t. 

Changing  from  s  to  t  in  Eq.  (2.5c),  (2.6),  and  (2.7) 
and  using  Eq.  (2.8),  we  see  that  these  equations  can  be  written 
in  the  form 


+  y(s(t)) sirij5(s(t))|  st(t)dt 


(2.9b) 


(2.10) 
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Equations  (2.8),  (2.9a  and  b)  ,  and  (2.10)  combined  with 
Eq.  (2.5a  and  b)  show  that  once  the  function  s(t)  has  been  deter¬ 
mined,  the  desired  pressure  distribution  coefficient,  Cp(s),  and 
the  drag, lift,  and  moment  coefficients  c^,  c^,  and  cm  can  be  in¬ 
ferred  by  direct  calculation. 

r  o  i 

We  now  invoke  one  of  the  main  results  in  Wu's  work1  , 
viz,,  that  the  unknown  function  s(t)  is  the  solution  to  a  non¬ 
linear  integral  equation,  viz.. 


J  (2-U) 

“1  I  “1  I  1 


where  f(t,tQ)  is  given  by  Eq.  (2.1),  with  the  real  constant,  A, 

and  the  complex  constant,  tQ,  still  to  be  determined.  It  will 

be  recalled  that  t  is  the  image  in  the  t-plane  of  the  point 

z-oo.  Inposition  of  the  condition  that  the  complex  velocity 
-ia 


-HJe 


as  z-*oo  leads  to  the  following  equation  for  tQ: 


t 

°  vis 


exp 


-1 


XT 


dr] 


(2.12) 


Finally,  if  S  is  the  (given)  normalized  length  of  the 
wetted  perimeter,  we  have 


s(l)  =  S  (2.13) 

The  function  s(t)  and  the  constants  A  and  CQ  are  to 
be  obtained  by  solving  Equations  (2.11),  (2.12),  and  (2.13) 
simultaneously.  It  is  this  solution  process  which  we  carry  out 
numerically  and  for  which  we  have  developed  the  computational 
procedure  described  in  the  next  section. 
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3.  COMPUTATIONAL  PROCEDURE 

To  calculate  the  solution  to  Equations (2. 11) -(2. 13) ,  we 
use  the  method  of  successive  iterations.  For  this  purpose,  it 
is  convenient  to  write  Equations  (2.11)  and  (2.13)  in  a  slightly 
altered  form.  We  put 

f(t,t0)  -  Ag(t,to)  (3.1) 

and 

s(t)  =  AN(t)  (3.2) 


where  g(t, tQ)  is  the  algebraic  function  multiplying  the  constant 
A  in  Eq.  (2.1). 

Equation  (2.11)  then  becomes 


while  for  Eq.  (2.13),  we  get 


(3.3) 


A 


S 

"  HU7 


(3.4) 


There  are  many  ways  in  which  to  start  off  the  iteration 

process.  We  have  followed  Wu  in  introducing  the  idea  of  the 

"basic  flow."  This  is  a  known  two-dimensional  cavitated  flow 

about  some  body  which  ideally  -  but  not  necessarily  -  has  a  close 

resemblance  to  the  given  body.  From  this  known  basic  flow,  one 

can  obtain  the  starting  quantities  or  zeroth  iterates,  s^(t) 

and  t  (°)*.  s^(t)  and  t  ^  are  substituted  into  the  right  side 

o  o 

of  Eq.  (3.3)  to  yield  N^(t),  from  which  is  then  calculated 

*Since  we  start  off  with  s^(t)  rather  than  N^(t),  we  do  not 
need  the  zeroth  iterate,  A^. 
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via  Eq.  (3.4).  Likewise  tQ^  is  obtained  by  substitution  of 
s°(t)  and  tQ(0)  into  the  right  side  of  Eq.  (2.12).  Finally, 
s(1)(t)  =  A^N^^(t),  in  accordance  with  Eq.  (3.2). 

The  formalities  of  going  from  the  (n-l)st  to  the  nfc^ 
iterate  are  defined  by  the  following  equations: 


(3.5) 


A<a>  =  -t-S  _ 
N^U) 

s(n>(t)  «  A^N^(t) 


(3.6) 

(3.7) 


(3.8) 


To  these,  we  may  add  equations  for  the  pressure  coef¬ 
ficient  and  the  drag,  lift,  and  moment  coefficients,  quantities 
which  although  not  iterated  in  themselves,  are  to  be  calculated 
at  the  end  of  each  iteration.  We  have  from  Eq.  (2.8) 


C 

CP 


(t) 


(3.9) 


while  from  Eq.  (2.5a  and  b),  (2.9a  and  b) ,  and  (2.10)  we  get 


cd(n)  *  cx^cosa  +  Cy^sina 

(3. 10a) 

■  *  cx^sina  +  Cy^cosa 

(3. 10b) 
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(3.11) 


X(8(n)(  t))C0s2(s(n)(t» 


+y(s(n)(t))sin3(s(nJ(t)) 


st<n)  (t)dt 


(3.13) 


As  of  this  writing,  we  have  not  Investigated  the  theore¬ 
tical  convergence  of  the  foregoing  sequences.  However,  such  con¬ 
vergence  seems  likely  when  the  choice  of  basic  flow  is  not  too 
unreasonable,  and  indeed  there  is  numerical  evidence  of  conver¬ 
gence  in  certain  of  the  specific  cases  to  be  described.  At  any 
rate,  we  shall  assume  that  if  the  zero**1  order  iterates  are  suit¬ 
ably  chosen,  the  sequences  defined  by  Equations  (3. 5) - (3. 12)  do, 
in  fact,  converge.  Note  that  even  when  this  assumption  is  cor¬ 
rect,  it  does  not  follow  that  the  numerical  realizations  of  these 
sequences  converge,  since  purely  numerical  phenomena  such  as  ac¬ 
cumulation  of  errors  may  intervene  to  spoil  the  convergence. 

The  actual  cct^mtation  divides  itself  logically  into 
four  stages: 

a)  numerical  specification  of  the  body  geometry, 

b)  basic  flow  calculation, 

c)  numerical  realization  of  the  iteration  procedure,  and 

d)  calculation  of  the  pressure  coefficient  and  the  drag. 
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lift,  and  moment  coefficients  from  Eq.  (3.7)  and  Eq.  (3. 9) -(3. 13) . 
We  describe  each  of  these  in  turn. 

A.  Specification  of  the  Body  Geometry 

The  body  geometry  enters  into  the  foregoing  equations 
through  the  functions  ]3(s)  ,  x(s),  and  y(s)  .  We  first  relate  these 
functions  analytically  to  the  usual  x-y  equations  of  the  body  pro¬ 
file.  We  then  define  a  set  of  data  points  along  the  perimeter  at 
which  F(s) ,  x(s) ,  and  y(s)  are  to  be  eval.uatcd. 

Let  the  body  shape  in  the  wetted  region  be  specified 


originally  by  two  single  valued  functions,  viz. , 

9  *  9+  $)  upper  profile  (9  >  0)  (3.14a) 

9=9"  (x)  lower  profile  (y  <  0)  (3. 14b) 

which,  upon  the  introduction  of  as  the  unit  of  length*,  become 

y  -  y+  (x)  y  >  0  (3. 15a) 

y  -  y"  (x)  y  <  o  (3.15b) 


These  functions  are  to  be  continuously  differentiable 
except  at  the  nose,  where  the  profile  has  a  vertical  tangent. 

Under  the  assumption  that  the  nose  is  parabolic,  it  is  convenient 
to  Introduce  a  new  variable,  w,  related  to  x  by  the  equations 

*  -  -  V5T  upper  profile  (y  >  0)  (3.16a) 

w  «  Y)T  lower  profile  (y  <  0)  (3.16b) 

We  then  have  instead  of  Equations  (3.15),  the  single 

equation 

y  *  h(w)  w,^  £  w  £  (3. 17) 

We  recall  that  D  «  6 €,  the  length  of  $-axis  intercepted  by  the 
body}  see  Figure  2.1. 
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where 


h(w)  =  y*(w2) ,  w  5  0  (3.18) 

The  quantities  wffl^n  and  w(nax  are  defined  by 

"min  "  '  V^T  <3- 19> 

"max  "  V*B  (3.20) 


where  x^  and  x^  are  the  x-coordinates  of  the  detachment  points 
A  and  B  shown  in  Figure  2.1.  In  these  equations,  the  positive 
square  root  is  intended. 

Because  the  nose  is  parabolic,  the  function  h(w)  is 
continuously  differentiable  throughout  the  interval  wm^n  <  w  < 
w  .  Thus  the  function  g(w) ,  defined  by 

*(»>  -  gjj,  (3.2i) 

is  continuous  in  that  same  interval. 

The  arc  length  8  can  be  expressed  as  a  function  of  w 
as  follows: 


•w  ■  /  V(M2  +  (*?- 

"min 


(3.22) 


Substituting  from  Equations  (3.16)  and  (3.21)  into  Equation  (3.22), 
we  get 


w 

M  -  J  V(2w’)Z  +  lg(w* ) )  2dw* 


min 


(3.23) 
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Finally,  from  the  definition  of  the  angle  £  of  Figure 

dx 

2.1,  it  follows  that  0  is  related  to  the  derivative  ^  by  the 


equations 


3  =  2v  -  arc  cos  g(w)  <  0,  w  ^  0 


_  „  dx 

3  =  arc  cos 


g(w)  >  0,  w  <  0 


3  =  2ir  +  arc  cos  g(w)  >  0,  w  >  0  (3. 

Here  the  arc  cosine  means  the  single  valued  inverse 
cosine  function  whose  range  extends  from  0  to  tr.  It  is  necessary 
always  to  express  3  in  terms  of  this  particular  branch  of  the  in¬ 
verse  cosine  since  it  is  only  this  branch  that  the  computer  sub¬ 
routine  calculates. 

From  Equations  (3.16)  and  (3.23)  we  substitute  for  dx 
and  ds  in  Equations  (3.24) -(3.26)  to  obtain 


(3.24) 


(3.25) 


(3.26) 


3  =  2ir  -  arc  cos 


r(2w;  +[g(w)  T 


,  g(w)  £  0,  w  £  0 


(3.27) 


3  -  arc  cos 


,  g(w)  >  0,  w  <  0 


(2w)2+[g(w)]‘ 


(3.28) 


3  =  2v  +  arc  cos 


,  g(w)  >  0,  w  >  0 


(2w)Z+[g(w)]- 


(3.29) 


By  means  of  Equations  (3.16),  (3.17),  (3.21/,  (3.23), 
and  (3. 27) -(3. 29)  we  have  expressed  each  of  the  variables  x,  y, 
s,  and  3  as  a  function  of  w  alone.  Thereby  x,  y,  and  3  are  de¬ 
termined  as  functions  of  s,  a  fact  which  we  use  later  in  the  pro¬ 
gram  when  carrying  out  the  iteration  procedure. 
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We  now  define,  along  the  wetted  perimeter  of  the  body, 
a  set  of  data  points  at  which  x,  y,  0,  and  s  are  to  be  evaluated 
Because  on  the  actual  profiles  of  interest,  0  changes  very  rapidly 
as  a  function  of  s  or  w  in  the  vicinity  of  the  nose  and  very  slowly 
elsewhere,  we  distribute  the  data  points  with  two  different  densi¬ 
ties.  In  the  vicinity  of  the  nose  a  high  density  or  fine  mesh  is 
used;  elsewhere  we  use  a  low  density  or  rough  mesh.  Let  and 
PpR  be  those  points  on  the  upper  and  lower  profiles  that  separate 
the  part  of  the  perimeter  over  which  0  changes  rapidly  from  the 
rest.  The  points  are  indicated  in  Figure  2.1,  where  their  ab¬ 
scissas  are  denoted  by  and  For  a  given  profile  these 

points  -  and  therefore  the  values  of  and  -  are  chosen  by 
eye.  In  the  discussion  which  follows  it  will  become  evident  that 
in  general  only  one  of  the  two  points  can  be  a  data  point  (mesh 
point)  ,  and  we  arbitrarily  take  this  one  to  be  Pjyp,  the  point  on 
the  upper  profile.  We  denote  by  PRR  the  mesh  point  lying  closest 

it 

to  PRR  on  the  side  toward  the  point  B.  The  quantitative  specifi¬ 
cation  of  PRR  will  be  given  presently. 

Corresponding  to  the  normalized  x-values  we  define  the 

w-velues 


WRF  "  '  VXRF  (3.30) 


Similarly,  we  can  write 


w 


FR 


(3.32) 


but  this  equation,  unlike  the  previous  two,  does  not  provide  a 
definition  of  the  quantity  on  its  left-hand  side  because  is 
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as  yet  unknown.  In  fact,  as  we  shall  see,  w?R  is  obtained  first, 

whence  can  be  calculated  from  Equation  (3.32),  if  desired. 

th 


w. 


min 


Let  w^  be  the  w-value  of  the  i  mesh  point,  counting 
as  w^.  Tnen  from  the  definitions  of  the  mesh  points 


and  we  have  the  scheme: 


wmin 

1  wi  <  WRf 

rough  mesh 

(3.  33) 

WRF 

1  wi  1  WFR 

fine  mesh 

(3.34) 

WFR 

^  wi  ^  wmax 

rough  mesh 

(3.35) 

Now  choose  a  value  for  NR,  the  number  of  rough-mesh 
intervals  on  the  upper  profile.  Since  «min  and  w^  are  already 
defined  by  Equations  (3.19)  and  (3.30),  the  choice  of  NR  fixes 
the  value  of  (Aw)R,  the  length  of  the  rough-mesh  interval  in  w. 
Specifically, 


(Aw) 


R 


WRF 

N 


-wmin 
T - 

R 


(3.36) 


With  (Aw)R  thus  determined,  we  see  that  NR,  the  number 
of  rough-mesh  intervals  on  the  lower  profile,  cannot  be  arbitrary 
but  follows  from  (3. 35)  and  the  definition  of  P^.  We  have  in 
fact 

* 

+  c,  0  £  r  <  1  (3.37) 

R 

which  says  that  NR  is  the  integer  part  of  the  left  side  of  Equa¬ 
tion  (3.37).  This  determination  of  NR  then  leads  at  once  to  the 
desired  expression  for  Wpg,  viz.. 


*FR  ”  "max  ’  NR^R 


(3.38) 
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Finally,  we  choose  a  value  for  Np,  the  "umber  of  fine- 
★ 

mesh  intervals  .  From  Equations  (3.^0)  and  (3.38)  we  can  then 
calculate  (Aw)_,  the  length  of  the  firc-mesh  interval.  We  have 


(Aw)f  =  — £ 


WFH'WRF 


(3.39) 


The  foregoing  equations  now  enable  us  to  write  the  fol¬ 
lowing  expressions  for  the  calculation  of  w^,  the  w-value  of  the 

i^  mesh  point.  We  recall  that  w*  =  w  .  =  -VxT. 

i  mm  ’  A 

wi  =  “  V^O-DCM*  i=l,2,...NR  (3.40) 

wi  =  WRF  +  <l  *  n£  ‘  l)(^w)F  (3.41) 

*  “V*A  +  NR^^R  +  (i-HR‘1)(Aw)i 


i  =  NR  +  1, 


<  +  nf 


»  wFR  +  (i  -  NR  -  Np  -  l)(Aw)R  (3.42 

*  -  +  NR(Aw)R  +  Nf(Aw)f 

+  (i  -  NR  -  Np  -  1)(Aw)r 

i  »  Mr  +  Np  +  1,  . Nr+Sf+Hr+1 

Next  we  evaluate  the  functions  x(w) .,  h(w) ,  s(w],  and 

]5(w)  at  the  points  w^  defined  by  Equations  (3.40) -(3.42) .  From 

2 

Equations  (3.16),  x(wj)  ■  w^,  while  h(Wj)  is  evaluated  via  Equa¬ 
tions  (3.17)  and  (3.18),  stwj  via  Equation  (3.23),  and  £(wj) 
via  Equations  (3.27) -(3.29) .  We  should  note  that  for  the 

When  the  body  is  symmetric  in  y  -  a  common  occurrence  -  it  is  con¬ 
venient  to  have  the  nose  point,  0,  be  one  of  the  amah  points, 
w^.  Accordingly,  we  always  make  Np  an  even  number. 
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profiles  in  whic1’  we  are  presently  interested,  the  functions 
jr^(x)  -  and  therefore  h(w)  -  are  simple  combinations  of  elemen¬ 
tary  functions  so  that  the  evaluation  of  h(w^)  presents  no  problem. 
The  function  g(w)  appearing  in  Equations  (3.23)  and  (3.27) -(3.29) 

JL 

and  defined  as  is  evaluated  by  first  differentiating  h(w) 
analytically  and  then  computing  the  resulting  expression  at  the 
desired  value  of  w. 

To  evaluate  the  integral  on  the  right  of  Equation  (3.23), 
we  first  observe  that 

i-1 

s[wi]  =  Y.  1  j  1  =  2’3»  ••  •  (3.43) 

J-l 

where  Ij  stands  for  the  integral  from  Wj  to  wj+i*  For  computa¬ 
tional  purposes,  it  is  advantageous  to  express  Equation  (3.43) 
in  recursive  form,  viz., 

s[w1+i]  =  slwjJ  +  1^  i  =  1,2 .  (3.44) 

with  s[w^3  =  0. 

For  any  given  i,  the  integral  I<  in  Equation  (3.44)  is 
calculated  repeatedly  by  means  of  Simpson's  rule,  the  number  of 
Siopson  subintervals  doubling  at  each  repetition.  The  process 
stops  when  successive  computed  values  of  differ  by  less  than 
some  preassigned  fractional  amount,  in  cur  case  0.1X. 

This  coapletes  the  description  of  Fart  A  of  our  pro¬ 
gram,  viz.,  the  specification  of  the  body  geometry  (profile  cal¬ 
culation).  We  conclude  with  a  recapitulation  of  the  input  and 
output  variables  for  Part  A. 


Inputs  for  Part  A 


t,  $A,  N+,  Np,  Equations  y  -  $+(x)  ,  y  =  y'(x) 


Outputs  from  Part  A 


wis  i  =  1,2,  ...  +  Np  +  +  1,  x(v^)  , 

s[wi],  P(w.). 

In  addition  to  being  stored  for  use  in  the  rest  of  the 
program,  these  five  variables  are  printed  out  to  expedite  program 
analysis  and  checking. 

B.  Basic  Flow  Calculation 

As  already  mentioned,  we  use  some  known  flow,  the  so- 
called  "basic  flow,"  to  furnish  the  zeroth  order  quantites,  s^°^(t) 
and  t  ^ ,  with  which  the  iteration  process  is  started.  The  choice 
of  the  basic  flow  is  made  on  the  dual  grounds  of  ease  of  calcula¬ 
tion  and  a  presumed  resemblance  to  the  flow  around  the  body  of 
interest.  In  this  subsection  we  exhibit,  by  way  of  illustration, 
one  such  basic  flow,  and  we  show  how  it  leads  to  the  numerical 
specification  of  the  zeroth  order  iterates. 

In  the  case  in  question  -  and,  for  that  matter,  in  all 
the  cases  thus  far  computed  -  the  original  body  is  a  parabola  lo¬ 
cated  symmetrically  with  respect  to  the  x-axis.  The  base  AB  of 
the  body  -  see  Figure  (2.1)  -  is  normal  to  the  x-axis  so  that 
XE  =  XA*  w^ile  frora  symmetry  we  have  yR  =-  y^.  We  assume  an 
incident  flow  parallel  to  the  5-axis,  i.e.,  a-0. 

For  a  basic  flow  we  choose  the  flow  about  a  wedge  with 
sides  of  length  2  located  symmetrically  with  respect  to  the  5-axis 
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and  having  its  vertex  at  0  (see  Figure  3.1).  This  flow  has  the 
merit  that  its  s(t)  function  is  known  explicitly  [2].  As  with  the 
parabola,  the  base  A'B*  of  the  redge  is  normal  to  the  'x-axis, 
while  the  incident  flow  is  parallel  to  that  axis.  The  fixed  de¬ 
tachment  ooints  of  the  basic  flow  are  at  A*  and  B’  corresponding 
to  the  values  t=-l  and  t=l,  respectively,  in  the  t-plane.  The 
half-angle  of  the  wedge  at  0  is  denoted  by  e-rr  with  0  <  €  <  y. 

To  coordinate  the  wedge  dimensions  and  the  basic  flow 
parameters  with  the  parabola  and  its  flow,  we  proceed  as  follows: 
we  set  the  length  of  the  base  A’B1  of  the  wedge  equal  to  that  of 
the  parabola,  AB.  This  gives  the  cavities  in  the  two  flows  the 
same  initial  widths.  Next  we  require  that  the  wetted  perimeters 
of  the  two  bodies  have  the  same  lengths,  i.e., 

-  (3.45) 


These  two  conditions  fix  the  angle  en  since  we  have 


err  =  arc 


sin 


(3.46) 


Further,  the  cavitation  number  a  is  assumed  to  be  the 
same  in  both  flows.  Then,  since  the  velocity  magnitudes  are 
normalized  to  the  free  streamline  speeds  in  both  flows,  the  nor¬ 
malized  incident  flow  speeds  must  be  equal  and  in  fact  are  given 
by 


U  » 


(3.47) 


Finally,  all  lengths  in  both  flows  are  normalized  with 
1  (=  0(5  of  Figure  2.1). 


respect 
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The  quantities  s^°^(t)  and  t  ^  ,  arising  from  the 
wedge  flow  and  serving  as  zeroth  order  iterates  in  the  parabola 
flow,  can  now  be  written  down  from  Wu’ s  solution  to  the  wedge 
problem  [2].  Specifically, 


(°) 


(t)  =  I  -  sw(-t) 


s(o)(t)  -  i  +  sw(t) 


-1  1  t  <  0 
0  1  t  £  1 


(3.48) 

(3.49) 


where  the  function  sw(t)  is  defined  by 


sw(t)=i 


t2(1"t)(t2+V2)  (t2+V~2)  +2 el  t' 1"2e(t'2+V2)  (t,2+V”2)1dt 


V2 (1+V2)  +2e j  t,1"2e(t'2+v2)1(t'2+v'2)1dt' 


(3.50) 


In  Equation  (3.50),  the  constant  V  is 

_  1 

V  =  (l+o)  **  (3.51) 


It  will  be  recalled  that  €,  the  coefficient  of  tt  in 
Equation  (3.46),  is  restricted  to  the  interval  0  <  €  <  j. 

Note  also  that  sw(t)  is  defined  only  for  t  ^0. 

From  Wu's  solution,  the  constant  t  ^  is  given  by 

tQ(o)  —  iV  (3.52) 

Next  we  represent  s^(t)  numerically.  For  this  pur¬ 
pose  we  establish,  over  the  interval  -1  £  t  £  1,  a  set  of  equally 
spaced  mesh  points  at  which  s^°^(t)  is  to  be  evaluated.  The 
number  of  such  points  is  chosen  equal  to  the  number  of  mesh 


24. 


points  in  w,  i.e.,  +  %  +  +  1.  Since  the  t-points  are 

equally  spaced  ,  it  follows  that  the  mesh  width  is  given  by 


(At)° 


2 

K  +  NF  +  NR 


(3.53) 


where  the  superscript  o  stands  for  zeroth  iteration. 

Denoting  the  i^  mesh  point  in  the  zeroth  iteration 
by  Ci<°>  with  t|^°^  =-  1,  we  have 

ti(o)  =  -  1  +  (i-l)(At)°  i=l ,2 ,  ...  n£+Nf+N"+1  (3.54) 

or,  in  view  of  Equation  (3.53), 

t.(o>  =  -  1  +  LlzlL.  i=l,2 ,  ...  N^+Np+Nl+1  (3.55) 

nr+nf+nr 

We  now  confute  s^0^(t^0^)  and  t  ^  via  Equations 
(3.48) -(3.50) ,  (3.52),  and  (3.55).  It  should  be  remarked  that 
the  evaluation  of  the  integrals  in  Equation  (3.50)  is  carried 
out  according  to  the  scheme  already  described  in  connection  with 
the  integral  of  Equation  (3.23). 

Inputs  for  Part  B 

L,  1,  N*,  Np,  N”,  o 


Outputs  from  Part  B 


t 


i 


i-1,2. 


. ..  N^+Np+Ng+1 , 


This  single -density  mesh  is  used  only  in  representing  the  zeroth 
order  iterate,  s(°/(t).  When  computing  s(n)(t),  n  >  1,  we  dis¬ 
tribute  the  mesh  points  over  the  interval  -1  £  t  £  I  with  two 
densities,  in  much  the  same  manner  and  for  the  same  reason  as 
was  done  with  the  w-points  of  the  profile  calculation.  This  is 
discussed  more  fully  in  the  next  subsection. 
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C.  Numerical  Realization  of  the  Iteration  Process 

We  have  already  mentioned  that  the  mesh  points  at  which 
the  n*'*1  iterate  s^(t),  n  ^  1,  is  to  be  evaluated  are  distributed 
on  the  interval  -1  £  t  <  1  with  two  densities  (rough  and  fine  mesh) 
in  analogy  with  the  mesh  points  in  w.  Unlike  the  w-points,  how¬ 
ever,  the  mesh  points  in  t  change  from  one  iteration  to  the  next; 
so  also  does  the  number  of  mesh  points.  Accordingly,  we  denote 
the  t-value  of  the  i^  mesh  point  in  the  nfc^  iteration  by  t^n^  , 

i=l,2,  ...  N,j/n)  (t^n^=-l)  ,  where  NT^  is  the  number  of  mesh  points 
t*h 

used  in  the  n  iteration. 

As  the  first  step  in  explaining  how  the  numbers,  t^n^  , 
are  arrived  at,  we  describe  the  procedure  for  calculating  t^p^ 
and  tpR^  ,  n  ^  1;  these  are  defined  as  the  mesh  points  in  t  at 
which,  with  increasing  t,  the  mesh  changes  from  rough  to  fine  and 
fine  to  rough,  respectively.  That  is,  they  are  the  analogies  of 
Wpj.  and  WpR  except  for  the  dependence  on  iteration  number.  The 
calculation  of  tpp^  and  tpR^n)  recurs^ve  since  it  presumes 
that  the  mesh  points  t^n~^  have  already  been  established. 

Let  s^p  and  SpR  be  the  arc  lengths  at  which  the  w-mesh 
density  changes.  We  may  write 

SRF  =  s  1  (3.56) 

8FR"s^Fr1  (3-57) 

where  it  will  be  recalled  that  the  values  of  the  function  on  the 
right  of  these  equations  form  part  of  the  output  of  Part  A.  It 
is  important  to  note  that  since  s(wp  depends  only  on  the  geometry 
of  the  body,  the  quantities  s^  and  sfR  are  independent  of  the 
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iteration  number.  Next  consider  the  (n-l)st  iterate  s^n“^(t). 

We  choose  as  t^P  the  value  of  t  that  satisfies  the  equation 

s(n_1)(t)  =  Srf  (3.58) 


Since  the  representation  of  s^n"^(t)  is  numerical  - 
it  is  defined  by  the  sequence  of  calculated  values  s^n“^  (t^n“^)  , 
i=l,2,  ...  NT^n”^  -  we  must  solve  Equation  (3.58)  by  interpolation 
among  these  calculated  values.  We  use  a  Lagrange  four-point  inter¬ 
polation  scheme  as  a  subroutine  for  the  numerical  solution  of 
Equation  (3.58). 


With  computed,  we  choose  a  value  for  the  number 

of  rough-mesh  intervals  in  t  on  the  upper  profile.  It  is  con¬ 
venient  to  make  this  number  equal  to  its  counterpart  in  the  case 

-f- 

of  the  w-mesh,  viz. ,  N^.  Note  that  this  particular  number  charac¬ 
terizing  the  t-rnesh  is  then  independent  of  n.  We  can  now  write 

t  (n)_t(n)  t  (n)+1 
(n)  _  CRF  tl  _  *KF 


(At), 


n: 


n: 


(3.59) 


where  (At)^  is  the  width  of  the  rough-mesh  interval  correspond¬ 
ing  to  the  nt^1  iteration. 

Before  calculating  tp^ ,  we  introduce  a  quantity 
tp^*.  This  is  defined  as  the  value  of  t  satisfying 

s  («-!)( t)  =  sFR  (3.59a) 

The  remarks  pertaining  to  Equation  (3.58)  apply  to  this 
equation  as  well.  In  particular,  the  solution  is  effected  via 
the  same  Lagrange  interpolation  process. 

Note  that  in  general  t^*  is  not  one  of  the  mesh 


points  in  t. 
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We  now  pick  the  mesh  point  tp^  as  follows:  starting 
at  the  point  t^  (=1)  ,  we  lay  off  in  the  direction  of  decreasing 
t  a  rough  mesh  of  width  (At)^,  where  (At)^  is  given  by  Equa¬ 
tion  (3.59).  From  the  mesh  points  thus  generated,  we  choose  as 
tp^  the  one  lying  closest  to  (If  tp^n^*  is  -  to  within 

the  number  of  significant  figures  retained  -  precisely  midway 
between  two  mesh  points,  then  the  mesh  point  with  the  less  posi¬ 
tive  t-  value  is  taken  to  be  tFRn^*)  If  N^n-  denotes  the  number 

th 

of  rough-mesh  intervals  on  the  lower  profile  in  the  n  itera¬ 
tion,  it  can  be  shown  that  the  analytic  representation  of  the 
procedure  just  described  is  given  by  the  equation 


t<n)  -  t^*  1  3  (n)* 

nt  FR  +  7  j  tpR  (n). 

- + 

0  1  r  <  1 


(3.60) 


where  the  first  equality  follows  from  the  fact  that  t^n^  =  1. 

The  second  equation  states  that  NRn^ ’  is  the  integer  part  of  the 
expression  on  the  left.  This  determination  of  N^”  leads  at 
once  to  the  desired  expression  for  tpR^  (cf.  Equation  (3.38) 
for  wFR) : 


t(n)  _  t(n)  <n>~rAt^n) 
CFR  "  cNt  -  Nr  'ac;R 


(3.61) 


=  1  -  N 


(n)- 
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Next,  we  choose  the  number  of  fine-mesh  intervals  in 

t  to  equal  N„,  the  number  of  such  intervals  in  w.  As  in  the  case 
r 

of  NR,  this  number  is  independent  of  n.  From  the  meaning  of  t^ 
and  t^  together  with  their  calculated  values  as  given  by  Equa¬ 
tions  (3.58)  and  (3.61),  we  can  then  write  for  (At)£n\  the  width 
of  the  fine -mesh  interval 

t(n)  -  t<n> 

(At)*n)  =  . —  (3.62) 

The  foregoing  equations  now  enable  us  to  use  the  fol¬ 
lowing  expressions  for  the  calculation  of  the  quantities  t^ . 

We  recall  that  t/n^  =  -1. 

t<n)  =  -1  +  (i  -  l)(At)*n)  i  =  1,2,  ...  (3.63) 

t<n)  =  t<£>  Mi-N+R-  l)(At)<n> 

=  -1  +  N^(At)^n)  +  (i  -  NR  -  l)(At)<n>  (3.64) 

i  -  Nj  +  1,  _  Nr  +  NF 

Cin)  =  Cre}  +  U  -  n£  -  Np  -  l)(At)<n) 

*  -1  +  NR(At)<n)  +  Np(At)£n) 

(3.65) 

+  (i  -  NR  -  Nf  -  l)(At)<n) 

1  "  NR  +  NF  *  l*  -  NTn) 

where  »  NR  +  Np  +  NR^n)”  +  1. 
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With  the  mesh  points  tj^  established,  we  can  discuss 
the  computation  of  the  nC^  order  quantities  s^(t^n^)  and  t^n^ 
via  the  iteration  Equations  (3.5)  -  (3.8).  The  principal  task 
is  to  perform  numerically  the  integrations  occurring  in  those 
equations.  This  is  a  complex  operation  having  a  number  of  aspects 
which  we  now  discuss. 

1.  Consider  first  the  integral  over  t' ,  the  "outer  inte¬ 
gral"  of  Equation  (3.5).  Since  the  function  (t)  -  and  there¬ 
fore  s<n>(t),  see  Equation  (3.7)  -  is  to  be  evaluated  at  the 
points  t|n^  given  by  Equations  (3.63)  -  (3.65),  it  follows  that 
the  upper  limit,  t,  of  the  outer  integral  is  to  be  assigned  the 
values  t^  only.  Now  it  fums  out  that  the  integrand  -  call  it 
I(n)(t’)-  of  the  outer  integral  is  a  smooth  function*  of  t’ ,  and 
therefore  we  save  a  great  deal  of  machine  time  without  much  loss 
of  accuracy  by  evaluating  I^(t')  only  at  the  points  t'  «  tj^ . 

The  integrals 

t(n) 

N(n)(t(n))  m  J  l<n>(t,)dt’  i  -  1,2,  ...  (3.66) 

-1 

are  then  computed  by  a  recursion  scheme  almost  identical  to  that 
defined  by  Equation  (3.44).  The  only  difference  is  that  in  the 


*  This  may  not  be  true  in  the  neighborhood  of  t'-O  for  the  10-to-l 
parabola  at  sero  angle-of-cttack.  See  the  discussion  of  Case  1 
in  Section  4. 
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present  case  the  integrand  is  given  in  tabular  form 

rather  than  analytically ;  thus  the  values  of  Iv  "(t’)  required 
by  the  Simpson's  rule  subroutine*  must  be  obtained  by  interpo¬ 
lation  among  the  quantities  l^(tj^).  For  this  purpose,  as 
with  all  our  interpolations,  the  Lagrange  four-point  procedure 
is  used. 

2.  We  can  write  I^(t’)  as 


I  ^  ( t* )  =  exp 


l-c'2  j(n>(t') 


7 r 


gt. 


(3. 67) 


where  J^(t')  represents  the  "inner  integral" 


X 

J(n)(t')  =  f  d T  (3.68) 


The  factor  gt.  (t,,tQ^n"1^)  p-  in  Equation  (3.67)  is  a 
simple  algebraic  function  of  t'  and  t^n~^  ,  as  can  be  seen  from 
Equations  (2.1)  and  (3.1);  its  evaluation  at  t'  «  therefore 


*  This  subroutine  has  already  been  alluded  to  in  the  discussion 
following  Equation  (3.44).  The  successive  refinement  of  the 
Simpson  mesh  should  stop  once. the.  density  of  the  o imp ton  points 
exceeds  that  of  the  points  tt^n"1'.  Further  refinement  pro¬ 
duces  no  improvement  in  the  approximation  to  the  integral 
since  the  additional  integrand  values  are  being  obtained  by 
interpolation  anyway. 
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presents  no  difficulty.  On  the  other  hand,  because  the  function 
$(s^n"^(r))  ordinarily  varies  very  rapidly  when  the  value  of 
S(n~*)(T)  corresponds  to  the  vicinity  of  the  nose,  a  special 
technique  hat  to  be  used  to  carry  out  the  -integration  in  Equa¬ 
tion  (3.68).  In  this  proc.dure  -  a  so-  ailed  adaptive  integra¬ 
tion  subroutine  -  Simpson' s  rule  is  applied  separately  to  vari¬ 
ous  subintervals  of  the  integration  interval  -1  £  t  <  1.  The 
Simpson  mesh-width  is  made  smallest  in  that  subinterval  where 
the  integrand  is  varying  most  rapidly,  etc.  The  details  of  this 
subroutine  can  be  found  in  [4]. 

3.  The  application  of  the  adaptive  integration  technique 
requires  that  certain  additional  interpolations  be  performed. 

The  integrand  function  ]3(s^n~^(T))  is  known  only  at  the  points 
r=.t^n" ^corresponding  to  the  values  (t£n~^)  outputted 

from  the  (n-l)st  iteration.  To  evaluate  F(s^n"^(x))  at  one 

of  the  Simpson  points,  we  proceed  as  follows:  denoting  by  t. 

J  I” 

any  one  of  the  Simpson  points  in  t,  we  first  compute  s^n”^(r.  ) 

J  I1* 

by  interpolation  among  the  known  quantities  s^n  ^(t^n~^). 

Denoting  the  resulting  value  of  s  by  sjn  ^ i  we  find  a  corre¬ 
sponding  value  of  w  through  the  equation 

s (w)  =  «<n-l>  (3.69) 

where  it  will  be  recalled  that  the  function  s(w)  depends  only  on 
the  geometry  of  the  body. 

Equation  (3.69)  is  solved  -  call  the  solution  wjn  ^  - 
by  interpolation  among  the  quantities  w^,  1  *  1,2,  ...  .  We  then 
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compute  the  function  fj (w)  at  w(n  ^  by  direct  evaluation  of  the 
formulas  given  in  Equations  (3.27)  -  (3.29).  The  resultant  num¬ 
ber,  P(w;n  ^ )  ,  is  the  desired  value  of  the  integrand  function  , 
0(s^n  ^ (t.  n)),  at  the  Simpson  point  r. 

Note  that  by  this  round-about  procedure  we  have  /'voided 
interpolation  on  a  rapidly  varying  function,  viz.,  f3(w)  .  [P(w) 

changes  rapidly  in  the  vicinity  of  the  nose,  w  =  0.  ] 

4.  Although  the  integrand  of  the  inner  integral  (see  Equa¬ 
tion  (3.68))  is  a  continuous  function  of  r  in  the  interval 

-1  <  t  <  1,  the  integrand  becomes  indeterminate  when  r  ~  t* . 
Evaluation  by  1' Hospital's  rule  proves  inconvenient  because 
the  necessary  derivatives  with  respect  to  r  cannot  be  computed 
with  sufficient  accuracy.  This  effect  is  particularly  marked 
when  t'  =  +1»  for  which  case  it  can  be  seen  that  the  second 
derivative  of  the  numerator  is  required.  Accordingly,  we  eval¬ 
uate  the  integrand  at  t  =  t1  by  assigning  to  it  the  value  at  a 
neighboring  point.  The  neighboring  point  is  chosen  close  enough 
for  the  integrand  to  differ  very  little  from  its  value  at  r  =  t1 
but  not  so  close  that  the  indeterminacy  spoils  the  accuracy  of 
the  evaluation.  Experimentally  it  has  been  found  that  a  separa¬ 
tion  of  about  0.003  between  t'  and  the  x-value  of  the  neighboring 
point  works  quite  well.  The  integrand  at  the  neighboring  point 
is  computed  by  interpolation  among  the  points  t![n”^. 

5.  The  integral  on  the  right  of  Equation  (3.8)  has  the 
same  form  as  the  inner  integral  J^n^(t’)  (Equation  (3.68))  except 
for  the  replacement  of  t' ,  which  is  real,  by  the  complex  quantity 
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t^n  ^ .  This  removes  the  aforementioned  indeterminacy;  otherwise 
the  same  considerations  apply  to  the  computation  of  this  integral 
as  in  the  case  of  J^(t'). 

We  conclude  this  subsection  by  remarking  that  the  itera 
tion  process  is  terminated  in  one  of  two  ways.  Either  a  maximum 
number  of  iterations  is  chosen  in  advance,  or  else  the  process 
stops  when  a  certain  convergence  criterion  is  satisfied.  The 
criterion  is  that 

t  \ 


Max 


Re 


fc<n)  -  t'-"'1)  1 


it(n)  „  t(n_1) 

o  o 


7W 


N 


(n)(i)-N(n-l)(1) 

N(n)(l) 


<  6  (3.70) 


where  the  number  6  is  generally  taken  to  be  10"^.  We  recall 
that  N^(l)  is  the  value  at  t  =  l  of  the  function  defined  in 
Equation  (3.5). 

Inputs  to  Part  C 


i  =  1,2, 


t(n-i) 

o 


n  =  1,2,  ...  (N<o)  =  +  NF  +  N"  +  1) 


Outputs  from  Part  C 

t(n),  S<n>(t<n>),  F(S(n)<t<n))),  i 

t  (n)  A(n)  n  =  1  2 

»  u  *■**■*••• 


=  1,2, 
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D.  Computation  of  Pressure,  Force,  and  Moment  Coefficients 
The  coeffi  . ients  C^n^[t],  ,  c^  ,  and  are 

computed  after  each  iteration  via  Equations  (3.9)  -  (3.13). 

can  effect  a  convenient  simplification  of  the  right  side 
of  Equation  (3.9)  by  differentiating  Equation  (3.5)  with  respect 
to  t  and  invoking  Equations  (3.1)  and  (3.7).  The  result  is  that 


u^r  =  texp 


.  1  -  tf  (n) 

7 r  J 


(t) 


(3.71) 


where  (t)  is  the  inner  integral,  Equation  (3.68). 
It  follows  that 


c<n)[t] 


2 

t  exp 


-7<l 


t2>j(n) 


(3.72) 


and  since  J^(t)  is  computed  at  the  points  t^  in  the  manner 
already  discussed.  Equation  (3.72)  immediately  furnishes  the, 


The  numerical  sped- 


quantities  i  =  1,2,  ... 

fication  of  the  function  C ^ (s)  -  the  normalized  pressure  as  a 
function  of  arc  length  -  then  results  from  associating  the  s-values, 
s<n>o4n>),  with  the  corresponding  C^r^[t£n^],  in  accordance  with 
the  identity 


C<n)(s(n)(t£n)))  =  C<n)[t[n)]  (3.73) 

The  integrands  in  the  integrals  of  Equations  (3.11)  - 
(3.13)  are  evidently  defined  numerically  at  the  points  t^.  We 
perform  the  integrations  by  combining,  as  previously  discussed, 


successive  applications  of  Simpson*  s  rule  with  interpolations 
among  the  integrand  values.  It  should  be  noted  that  the  con¬ 
tents  of  the  curly  brackets  in  the  integrands  are,  for  evalua¬ 
tion  purposes,  replaced  by  the  right  side  of  Equation  (3.72), 
while  s<n> (t)  itself  is  evaluated  by  solving  for  it  in  Equa¬ 
tion  (3.71).  Finally,  the  functions  x(s^(t))  and  y(s^  (t)) 
appearing  in  Equation  (3.13)  are  computed  by  interpolations 
among  the  tables  of  x(wi) ,  ,  and  s(w^)  that  form  part  of 

the  output  of  Part  A. 

Inputs  to  Part  D 


a,o 


), 


i  =  1,2, 


n  =  1,2,  . . . . 


Outputs  from  Part  D 


C<n>Lt<n>], 


i  =  1,2, 


n  =  1,2, 


•  •  #  • 
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4.  NUMERICAL  RESULTS 

In  the  Introduction,  we  summarized  the  results  obtained 
by  applying  the  foregoing  compu rational  scheme  to  four  ca vita ted 
flows.  We  now  present  some  of  the  numerical  details.  All  compu¬ 
tations  were  performed  on  the  IBM  7094  computer. 

Case  1.  10-to-l  parabola,  zero  angle  of  attack,  zero  cavita¬ 

tion  number. 

Inputs:  a  =  0,  a  =  0,  xA  =  xfi=  1,  =  xfTR  =  0.015, 

n£  =  30 ,  NF  =  100 ,  y^x)  =  ±  VT, 

Maximum  number  of  iterations  =  4.  Basic  flow=flow  about  wedge. 

At  the  present  time,  this  case  is  the  only  one  for  which 
there  seems  to  be  no  possibility  of  imputing  any  physical  signifi¬ 
cance  to  the  sequence  of  computed  pressure  distributions.  As  may 
be  seen  by  referring  to  Figure  4.1,  the  distribution  corresponding 
to  the  first  iteration  is,  to  a  high  degree  of  accuracy,  symmetric 
with  respect  to  the  nose  .  This  is  as  it  should  be  since  both  the 
basic  flow  and  the  true  pressure  distribution  are  symmetric.  On 
the  other  hand,  the  second  and  third  distributions  increasingly 
show  a  loss  of  symmetry,  while  the  fourth  distribution  bears  no 
resemblance  to  the  first  three  and  in  fact  is  almost  antisymmetric 
with  respect  to  the  nose.  The  corresponding  drag  coefficients, 
shown  in  Table  4.1,  vary  from  one  iteration  to  the  next  in  what 
seems  to  be  a  random  manner.  (The  lift  and  moment  coefficients 
vanish  when  the  angle-of-attack  is  zero.) 

The  nose  in  this  case  is  located  at  s=1.003  since  the  total  normal¬ 
ized  arc  length  for  a  10-to-l  parabola  equals  2.006.  Note  that  the 
interval  of  s -values  in  Figure  4.1  corresponds  to  the  vicinity  of 
the  nose,  where  the  very  rapid  variations  in  pressure  take  place. 
Outside  this  interval  C  goes  monotonically  to  zero  in  all  four 
iterations,  as  s  +  0  or  2.006. 


1  .DO 
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Ok.  * 


a  aw 


2nd  ITERATION 


Ist  ITERATION - 


4*n  ITERATION 


S*  ARC  L 


.96  .98  1.00 

GTH  FROM  DETACHMENT  POINT  A 


37. 


TABLE  4.1 


Iteration  number,  n 

<d(n) 

1 

0.00421 

2 

0.00483 

3 

0.00301 

4 

0.00409 

A  possible  explanation  of  these  results  is  being  sought 
in  the  fact  that  for  a  10-to-l  parabola,  the  integrand  of  the  inner 
integral  varies  extremely  rapidly  with  r  in  the  neighborhood  of  the 
nose,  r=0,  owing  to  the  behav  .or  of  H(s(t))  in  that  neighborhood. 
The  behavior  is  such  that  the  major  Lange  in  p  takes  place  over  a 
r-interval  of  width  smaller  than  that  o  ;  the  fine  mesh  interval. 
This  abrupt  variation  could  have  two  consequences:  first,  it  might 
induce  errors  in  the  evaluation  of  the  inner  integral;  second,  it 
can  be  shown  that  the  exact  value  of  the  inner  integral,  regarded 
as  a  function  of  t' ,  peaks  sharply  when  t'  is  in  the  region  where 
p  varies  rapidly.  Such  peaking  could  then  lead  to  inaccuracies  in 
the  computation  of  the  outer  integral,  owing  to  the  interpolation 
used  to  evaluate  the  outer  integrand,  I^(t'). 

If  must  be  emphasized  that  these  considerations  are  hypo¬ 
thetical.  The  effectiveness  of  our  computational  scheme  when  ap¬ 
plied  to  thin  bodies  having  rounded  noses  of  very  high  curvature 
is  subject  to  further  investigation. 

Case  2.  1-to-l  parabola,  zero  angle-of-attack,  zero  cavita¬ 
tion  number. 

Inputs:  a  =  0,  o=0,  xA  =  xB=l,  =  xpR  *  0.04, 

NR  =  15 ,  Np  =  40 ,  y~(x)  =  +  £  V** 

Maximum  number  of  iterations  *  7.  Basic  flow  *  flow  about  wedge. 
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Results  for  this  case  are  presented  in  Figures  4.2  and 
4.3,  and  Table  4.2.  Figure  4.2  shows  the  computed  pressure  coef¬ 
ficient  Cp  as  a  function  of  arc  length  along  the  parabola  for  each 
of  the  seven  iterations.  Since  the  flow  is  symmetric  about  the 
nose,  and  the  iterates  reflect  this  fact,  it  suffices  to  give  the 
pressure  distribution  over  half  the  perimeter.  Note  that  arc 
length  in  this  case  is  measured  from  the  nose  rather  than  from  the 
upper  detachment  point.  The  convergence  of  the  iteration  process 
is  suggested  by  the  fact  that  the  sixth  and  seventh  pressure  dis¬ 
tributions  differ  by  at  most  a  few  percent  (on  the  graph  they  are 
not  even  distinguishable) . 

In  Figure  4.3  we  show  a  comparison  between  our  final 
(seventh)  distribution  and  the  distribution  computed  from  an  ap¬ 
proximate  formula  derived  by  Johnson^.  Johnson's  formula  is 
meant  to  apply  to  parabolas  with  a  higher  fineness  ratio  than  one- 
to-one,  but  the  agreement  is  still  quite  close.  For  the  abscissa 
in  Figure  (4.3)  we  have  chosen  y  rather  than  s  because  the  area 
under  the  resultant  curve  can  be  shown  to  equal  half  the  drag  co¬ 
efficient  c^;  this  expedites  the  comparison  between  our  drag  coef¬ 
ficient  and  Johnson's,  which  turn  out  to  equal  0.307  and  0.289  re¬ 
spectively,  a  difference  of  less  than  6  percent. 

Table  4.2  shows  how  the  computed  value  of  c^n^  varies 
with  the  iteration  number,  n.  The  sixth  and  seventh  values  differ 
by  less  than  0.2  percent. 
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TABLK  4.2 


Iteration  No. ,  n 

c  <n> 
u 

1 

0.44199 

2 

0.31008 

3 

0. 30353 

4 

0.30950 

5 

0. 30597 

6 

0.30734 

7 

0.30673 

Case  3.  1-to-l  parabola,  1°  angle-of-attack,  zero  cavitation 
number. 

Inputs:  a=l°,  o  =  0,  xA=xB=lf  x^  =  0.25,  xpR  =  0.10, 

Nr  =  15,  Nf  =  120,  y^“(x)  =  + 

Maximum  number  of  iterations  =  14.  Basic  flow  =  flow  corresponding 
to  zero  angle  of  attack  on  same  parabola. 

This  case  is  distinguished  by  the  fact  that  for  the  func¬ 
tion  s^°^(t)  with  which  the  iterations  are  started,  we  choose  the 
(71 

final  iterate,  sv  1  (t)  ,  obtained  in  the  zero  angle-of-attack  case. 
Figure  4.4  shows  the  resulting  pressure  distributions  corresponding 
to  various  iteration  numbers  as  a  function  of  arc  length  measured 
from  the  upper  detachment  point.  Note  that  in  the  case  of  a  1-to- 
1  parabola,  S  =  2.323.  For  the  sake  of  clarity,  only  the  first, 
sixth,  thirteenth,  and  fourteenth  distributions  are  shown.  The 
still  considerable  separation  between  the  thirteenth  and  four¬ 
teenth  distributions  indicates  that  the  rate  of  convergence  is 
much  slower  than  in  the  zero  angle-of-attack  case.  It  can  be  seen- 
the  curves  not  shown  bear  this  out  -  that  the  pressure  distributions 
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shift  back  and  forth  as  the  iteration  number  changes.  A  manifes¬ 
tation  of  this  effect  is  exhibited  in  Figure  4.5  where  we  have 
plotted  the  s -value  of  the  stagnation  point  (maximum  pressure)  as 
a  function  of  iteration  number.  The  oscillatory  behavior  is  quite 
apparent,  as  is  the  gradual  approach  to  a  limiting  position. 

Table  4.3  gives  the  computed  values  of  ,  c^n^  , 

and  c  ^  as  a  function  n. 


TABLE  4.3 


Iteration  No. ,  n 

c  <»> 
d 

c  (n) 
ci 

c  <»> 
m 

1 

0.30727 

-0.00682 

-0.00042 

2 

0.30516 

0.03830 

0.01172 

3 

0.30298 

0.19697 

0.10083 

4 

0. 30795 

-0.01678 

0.00214 

5 

0 . 30488 

0.05413 

0.02038 

6 

0.30513 

0. 15519 

0.08022 

7 

0.30739 

0.01061 

0.01005 

8 

0.30577 

0.06780 

0.02979 

9 

0.30573 

0.12334 

0.06427 

10 

0.30786 

0,02816 

0.01722 

11 

0.30571 

0.07346 

0.03430 

12 

0.30965 

0.10423 

0.05441 

13 

0.30724 

0.04170 

0.02297 

14 

i 

0.30585  ! 

0.07581 

0.03649 

The  drag  coefficient  remains  substantially  constant, 
apart  from  a  small  random  fluctuation  -  rms  relative  amplitude 
about  1  percent  -  apparently  due  to  numerical  errors.  The 
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average  value  of  is  about  the  same  as  the  final  value  in  the 

zero-degree  case.  The  more  severe  fluctuations  in  c^n^  and  cm^ 
can  be  explained  by  referring  to  Equations  (2. 5a) -(2. 6b) .  From 
those  equations  we  see  first  that  when  a  is  small,  c^fc*  cx  and 
c I***  Cy.  Thus  the  computability  of  and  c^  depends  on  the  com¬ 
putability  of  the  integrals  in  Equations  (2.6a)  and  (2.6b)  respec¬ 
tively.  Observing  that  as  the  integration  variable  s  runs  from  0 
to  S,  the  angle  ]5(s)  goes  from  slightly  greater  than  tt  to  slightly 
less  than  2ir ,  we  see  that  sin^(s)  remains  negative  over  the  inter- 
val  0  <  s  S.  Since  Cp(s)  is  positive  on  that  interval  ,  the 
integrand  in  Equation  (2.6a)  is  negative  throughout;  thus  there  is 

no  cancellation  of  one  part  of  that  integral  by  the  other.  On  the 

3tt 

other  hand,  the  function  cos0(s)  is  antisymmetric  about  3  =  -£■ 

(s  =  ^S)  while  Cp(s)  is  almost  symmetric  about  the  same  point. 

It  follows  that  the  integral  in  Equation  (2.6b)  is  the  difference 
between  two  almost  equal  quantities,  and  its  computation  is  there¬ 
fore  highly  susceptible  to  numerical  error.  This  explains  why  the 


fluctuations  in  c 


(n) 


are  more  severe  than  in  c. 


(n) 


The  same 


argument  applies  to  c  ^  upon  remarking  that  the  terms  x(s)cosfJ(s) 
and  y(8)8ih^(s)  in  the  integrand  of  Equation  (2.5c)  are  both  anti¬ 
symmetric  about  s  =  -jS. 

It  should  be  noted,  however,  that  despite  the  foregoing 
effect,  some  tendency  toward  convergence  may  be  discerned  in  the 


tabulated  values  of  c.^  and  c  ^  .  In  fact,  the  numbers  seem  to 

m 


indicate  that  as  a  rough  guess,  the  limiting  value  of  c^  lies 


*For  a  few  values  of  n,Cp^(s)  goes  slightly  negative  over  a  short 
interval  in  s,  but  this  does  not  affect  the  argument. 
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somewhere  between  0.03  and  0.08,  while  for  c  it  lies  between  0.03 

tn 

and  0.06.  The  need  for  greater  accuracy  in  computing  these  inte¬ 
grals  is  evident. 

Case  4.  i.-to-l  parabola,  5°  angle-of-attack,  zero  cavitation 
number. 

Inputs:  a=5°,o=0,  xA=xR=l,  x^  =  0.25,  xpR  -  0.10, 

NR  =  15,  Np  =  120,  y^Cx)  =  + 

Maximum  number  of  iterations  ••=  14.  Basic  flow  =  flow  about  wedge. 

In  Figure  (4.6)  we  have  plotted  the  computed  pressure 
distributions  corresponding  to  the  first,  fourth,  ninth,  thirteenth, 
and  fourteenth  iterations,  the  remaining  distributions  having  been 
omitted  for  the  sake  of  clarity,  as  in  the  1°  case.  By  ccntrast 
with  that  case,  however,  there  seems  to  be  no  trend  toward  a  limit¬ 
ing  distribution,  and  our  present  opinion  is  that  the  sequence  of 
pressure  distributions  is  not  convergent.  On  the  other  hand,  we 
do  not  feel  that  these  distributions  are  meaningless  since  with  the 
exception  of  the  third  distribution  (not  shown) ,  which  has  a  large 
negative  dip,  they  are  all  of  the  expected  form  (judging  from  the 
0°  and  1“  cases),  while  the  computed  pressures  fall  within  the  ex¬ 
pected  range.  We  conclude  that  the  failure  to  converge  is  probably 
the  result  of  fluctuations  arising  from  numerical  inaccuracies.  As 
mentioned  in  the  Introduction,  this  idea  is  supported  by  the  com¬ 
parative  constancy  -  rms  fluctuation  over  fourteen  iterations 
about  5  percent-of  the  drag  coefficient  as  a  function  of  iteration 
number,  n,  and  by  the  quite  reasonable  figure  of  0.299  for  the 

- 

These  were  chosen  as  representative  of  all  fourteen  distributions 
in  the  range  of  pressures  encompassed. 
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average  over  the  iterations.  The  tabulation  of  these  values  and 
the  corresponding  values  of  the  lift  and  moment  coefficients  is 
exhibited  in  Table  4.4. 


TABLE  4.4 


Iteration  No. ,  n 

c 

cd 

c  (n) 
ci 

c 

m 

1 

0.30568 

>0.02721 

-0.00028 

2 

G. 28080 

0.14818 

0.05413 

3 

0.22256 

1.04143 

. 

0.50903 

4 

0.32191 

!  0.17081 

0.10733 

5 

o.  9101 

0.36652 

0.17795 

o. 2*0:4 

0.34706 

0.17524 

7 

0.29248 

0.50943 

0.25906 

8 

G. 30700 

0.22710 

0.11964 

9 

0.28604 

0.45990 

0.22614 

10 

0.31108 

0.29812 

0.16118 

11 

0.29963 

0.27486 

0.13723 

1? 

0.28485 

0.42457 

0.20862 

13 

0.30755 

0.39021 

0.20712 

14 

0.30790 

0.19595 

0.10319 

Note  the  much  larger  values  of  c^/n^  and  cm^  in 
the  case  of  the  5°  angle-of-attack  as  compared  with  the  1°  case. 
This  was  to  be  expected  and  further  substantiates  our  opinion  that 
the  5r  pressure  distributions  have  some  physical  significance. 
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5 .  SUPPLEMENTARY  DISCUSSION  AND  RECOMMENDATIONS 

As  evidenced  in  this  report,  our  experience  with  the 
computation  of  pressure  distributions  by  application  of  Wu' s 
theory  of  cavitated  flow  suggests  several  steps  that  might  be 
taken  to  improve  our  technique. 

1)  Since,  in  our  opinion,  the  numerical  inaccuracies  in 
the  various  integrals  are  probably  the  principle  cause  of  those 
anomalous  results  already  described,  a  major  effort  should  be 
made  to  improve  the  accuracy  of  the  i.nr.egration  process.  The 
most  direct  way  of  doing  this  is  simply  to  increase  the  number 
of  data  points,  especially  in  the  fine  mesh  region.  The  only 
limitation  on  the  improvement  in  accuracy  thus  attainable  is  the 
possibility  of  excessive  computation  time.  However,  tu  •  cases 
already  run  have  used  less  than  0.01  hours  per  iter  .ion  on  the 
IBM  7094,  so  that  we  seem  to  have  considerable  leeway  in  this 
respect. 

2)  It  is  possible  that  in  the  5°  case  there  is  an  addi¬ 
tional  factor  tending  to  prevent  convergence.  Specifically,  the 
basic  flow,  which  in  that  cas*.  was  the  flow  about  the  wedge  at 
0°  angle-of-attack,  may  be  too  far  removed  from  the  actual  flow. 
This  suggests  that  one  might  proceed  stepwise  by  using  the  0° 
solution  for  the  parabola  as  the  basic  flow  in  the  1°  case,  the 
1°  solution  as  the  basic  flow  in  the  2°  case,  etc.,  thereby 
"creeping  up"  on  the  situation  of  a  5°  (or  higher)  angle-of- 
attack.  Naturally  the  increments  in  the  angle-of-attack  could  be 
made  smaller  if  this  was  thought  desirable. 
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3.  An  increase  in  the  rate  of  convergence  of  the  iteration 
process  might  be  effected  through  some  technique  for  processing  the 
iterates  in  a  manner  more  sophisticated  than  successive  substitution 
into  the  integral  of  the  integral  equation.  For  example,  a  simple 
and  familiar  technique  of  this  kind  consists  of  substituting  into 
the  integral  not  the  last  iterate,  but  the  average  of  the  last  two. 
Apparently  this  has  been  found  to  work  quite  well  in  some  cases. 
Variations  of  this  idea,  as  well  as  more  elaborate  procedures,  could 
also  be  tried. 
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